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ABSTRACT 

High accuracy numerical quadrature methods for integrals of singular 
periodic functions are proposed. These methods are based on the appropriate 
Euler-Maclaurin expansions of trapezoidal rule approximations and their 
extrapolations. They are used to obtain accurate quadrature methods for the 
solution of singular and weakly singular Fredholm integral equations. Such 
periodic equations are used in the solution of planar elliptic boundary value 
problems, elasticity, potential theory, conformal mapping, boundary element 
methods, free surface flows, etc. The use of the quadrature methods is 
demonstrated with numerical examples. 
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1. INTRODUCTION 


In this work we shall present some Romberg-type quadrature methods for the 
numerical solution of singular and weakly singular Fredholm integral equations of 
the first and second kinds 

b 

of (t) + f K(t,x)f {x)±E = g{t), a^t<,b , (1.1) 

a 

where cj= 0 and a=l for first and second kinds respectively, and the corresponding 
eigenvalue problems 

b 

f (t) + \f K(t ,x)f {x)dx = 0, a <t<b. (1.2) 

a 

that arise, for example, from some boundary value problems over two dimensional 
domains with smooth and rectifiable boundaries. These quadrature methods are 
ultimately based on appropriate Euler-Maclaurin expansions for the trapezoidal 
■ rule, and treat the singularity in the kernel K(t,x) systematically. As such, these 
methods are simple, easy to implement, and have high order of accuracy. 

We assume that K(t,x) is periodic in both t and x with period T=b—a, and 
g(t) is periodic in t with period T. We furthermore assume that K(t,x) is of the 
form 


M 


K(t ,x ) = £ .*)i*-* r*(i°gu-z i ) ?fc 


Jb = l 




(1.3) 


where a k are real numbers satisfying a k > — 1, and p k are non-negative integers, 
H k (t t x), 1 < A; < M t and Hj(t ,x), j =1,2, are differentiable in t and x as many times 
as needed. In (1.3) it is assumed that H k (t t t ) ^ 0 whenever H k (t ,x) for some k. 
It is also assumed that ^ 0 whenever /^(f.x) ^0, and in this case the 

integrals in (l.l) and (1.2) are to be taken as Cauchy principal value integrals. 
When H k (t ,x) = 0, 1 < A: < M, and H x (t ,x) ^0, the integral equations (l.l) and (1.2) 
are called singular, and when H k {t t x) ^0 for some k 9 1< k < M f and H x (t t x) = 0, 
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they are called weakly singular. 

It is worth mentioning that the important cases of K(t,x) that arise in appli- 
cations are K(t ,x) = H i(t ,x)log \t-x | + Hz{t,x), K{t,x) = — 73 — — +H z (t,x), and 

b X 

combinations of the two. 

We shall assume that /(f), as K(t ,x) and g(t), is periodic in f , with period T , 
and is as many times differentiable as needed. We note that when (l.l) and (1.2) 
arise from some boundary value problems in two dimensions, defined over 
bounded domains with smooth boundaries, K(t,x), g(t), and /(f) also, in general, 
turn out to be smooth. 

One of the methods for solving (l.l) and (1.2) numerically is the quadrature 

method (see Baker (1977, Chap. 4, Section 3)), in which one replaces the integral 
6 

J K(t ,x)f (x)dx by a numerical quadrature formula, whose abscissas are 
' a 

1 n, with t -x it i=l,...,7L, then replaces the f (xj) by their approximation 

rv* pj 

fj, and finally solves the resulting system of linear equation for the fj. Obviously 
the accuracy of this method depends on the accuracy of the numerical quadra- 
ture formula being used, which in turn depends on the analytic properties of both 
the kernel K{t 9 x) and the solution f (t) over [a, 6 ]. It can be said, in general, that 
whenever K(t w x) is weakly singular or singular, the solution /(f) will be singular 
at the end points a and 6 . The singularity structure of /(f) may be complicated 
and difficult to determine; see MacCamy (1958) and Graham (1982) for some gen- 
eral results on this problem. When K{t ,x), g{t ), and / (f ) are (periodic) as 
assumed in the present work, then a and b in (l.l) and ( 1 . 2 ) can be replaced by a* 
and b ' respectively, where 6 '— a* = T. If we now assume that /(f) has singularities 
at a and 6 , then it should be singular at a* and 6 * ad hence at all f. As a result 
we conclude, heuristically, that /(f) cannot have any singularities, and this is the 
assumption that we have made above. 
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Let 

Xj = a + jh t h = (6-a)/n, n a positive integer. (1.4) 

Using the Euler-MacLaurin expansion for smooth integrands, and their extension 

to integrands having end point singularities (see Navot (1961, 1962)), in the next 

b 

section we derive Euler-MacLaurin expansions for the integrals f K(t,x)f (x)dx t 

a 

with K(t,x) as given in (1.3), We basically derive asymptotic expansions, for h-* 0, 
for the differences 

k{t,h) = I[t;f ] - I n [t;f ] f (1.5) 

where 

b 

/[*;/] = f K{t,x)f (x)di: (1.6) 

a 

and 

= s .*,-)/(*;). (1.7) 

j = i 

such that t is one of the points Xj, and is being held fixed, and uu n (f ,Xj) = hK{t ,Xj) 
for Xj ^ t, and depend on the type of singularity that K(t,x) has for t = x m 

Using the asymptotic expansions for A(£,/i), in Section 3 we derive Romberg-type 
numerical quadrature formulas, thus increasing the accuracy of /*[£;/] by as 
many orders of magnitude as we wish. In Section 4 we present quadrature 
methods that are based on these Bomberg-type formulas. In Section 5 we illus- 
trate the efficiency of our quadrature methods with numerical examples. In Sec- 
tion 6 we review some quadrature methods that have been proposed and bear 
some relation to the ones proposed in the present work. 

Note : The treatments of singular and weakly singular integral equations are 
not separated from each other throughout this paper. The reader interested in 
the treatment of singular integral equations could follow it easily by going through 
Theorems 2.1, 2.4, 2.7a, 3.1, 3.2, and the first part of Section 4, without having to 
consider the rest of this work. 
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2. EULEK-MACXAURIN EXPANSIONS 

The notation described below will be used throughout the remainder of this 
work. 

Let Xj = a+jh t j ^.l,...,™, h-(b —a)/n, where n is a positive integer. Let 
t e (a,6) be fixed and t e \xj | l<j<n— 1} for some n=n 0 . Obviously, there exists 
an infinite sequence of integers \n k ] k= 0 , 7i fc+1 >n k , A: =0,1,..., such that t is one of 

the Xj whenever n=n k , k =0,1 With the exception of Theorems 2. 1-2.3, the 

notation A-»0 will be assumed to mean that n->°° through the sequence of integers 

mg m 2 

We shall take 2' a j ( or S* a i ) to mean that a TOg (or a mi ) is to be multi- 

i=m t j=m l 

m z 

plied by 1/2, while a j vill he taken to mean that both a mi and are to be 

j=m 1 

multiplied by 1/2, in the respective summations. 

Theorems 2. 1-2.3, which are stated below without proof form the basis of our 
development throughout the remainder of this work. 


Theorem 2.1: Let the function g(z) be 2 m times differentiable on [a # 6]. Then 


D{h) = fg(x)dx - ^"g (xj) 

a j =0 


= (f^)! ~[ g(2M °( °(b)]^ 2M + R Zm jgf;(a,6) , 


where 


B, m , gW ., (x)fc (2.2) 

Here B M are the Bernoulli numbers, and B^x) is the periodic Bernoullian function 
of order fj . ;. In addition, since B^(x) are bounded on (-‘ 00 , 00 ), it follows that 

l^ 2 m[ff:(“. 6 )]l ^ M Zm (b-a)h 2m ma^\g^ m) {x)\, (2.3) 

where 
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M Zjn = max | B Zm (x) - B 2m | / (2m)! , (2.4) 

and, therefore, is independent of h. Consequently, if g(x) is infinitely 
differentiable on [ot,b], then D{h) has an asymptotic expansion of the form 

D{h)~ § (2m_1) (“) “ 9 l2M_1) ( b )]^ 2m as h -*0 . (2.5) 

For a proof of this result see Steffensen (1950). The expansion in (2.5) is the 
classical Euler-Maclaurin expansion for trapezoidal rule approximations of 
integrals of smooth functions. Navot (1961, 1962) has extended the Euler- 
Maclaurin expansions to trapezoidal rule approximations of integrals of functions 
having algebraic and/or logarithmic end point singularities. By a different 
approach that utilizes generalized functions, Lyness and Ninham (1967) have 
rederived Navot* s results. The results stated as Theorems 2.2 and 2.3 below, are 
special cases of those proved by Navot. 


Theorem 2.2: Let g{x) be 2771 times differentiable on [a,6] and let 

G(x) = (x-a) s sr(x), s>- 1. Then 

r A 

D(h) = J G(x)dx -/i2'G(xj) 

a j=i 

= ~%^r cfu ~' Xb)h ^ <a6) 

^=0 r* u 


where £(r) is the Riemann zeta function initially defined for Re r > 1 by 

oo 

£(r) = Yj and then continued analytically, and 
fc = l 


Pzm 



as 0 . 


(2.7) 


Itg(x) is infinitely differentiable on [a,b], then D(h) has an asymptotic expansion 


of the form 
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W~-S7 fSr G(2M_1) ( & )^ 

P= J ( T' * (2.8) 

- § A 5 W g M( a )/ t /i+»+i as h_o. 
m=o M- 

Starting from Theorem 2.2, Navot (1962) shows that the extensions of the 
Euler-Maclaurin expansion to trapezoidal rule approximations of integrals of func- 
tion of the form (a:-a) s [log(^:-a)] p flr(a:), with p being a positive integer and g(x) 
being sufficiently smooth in [ 0 , 6 ], can be obtained by differentiating both sides of 
( 2 . 6 ) p times with respect to s. Forp=l the following results are obtained. 


Theorem 2.3: Let g(x) be 2m times differentiable on [a. 6 ], and let 

G(x) = (x-a)*log(a? -a )£(£), s>-l. Then 


D(h) = f G(x)dx - h^'G(xj) 

a j=l 

= (2 - 9) 

- £ [-f(-S-M) + f(-s-#»)logA g0 *^ a ) A **** 1 +fs m , 

M =o L 

where f*(r) = d${r)/ dr, and 

~Pz m = 0^i 2m ) as h-*0. (2.10) 

If g(*) is infinitely differentiable on [a, 6 ] t then D(h) has an asymptotic expansion 
of the form 


z?(A) ~“? 1 (^r GC2M ' 1){6)A2/1 


M=1 

- £ |—C(— -*—/*) + f(-s-p,)iog h 


g^Ha) 


(2.11) 


f 1 - 


■h^ 3 


+ 1 


Note that both Theorem 2.2 and Theorem 2.3 are true for any s>— 1 , although 
they were originally stated for -l<s^ 0 . 

6 

x—t 


We shall now apply Theorems 2 . 1-2.3 to integrals of the types and 
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6 

f | x-t | s (log|x — 1 1 y 3 g (x)cLc , the former being defined as a Cauchy principal value 

a 


integral. 


Theorem 2.4: Let g{x) be 2m times differentiable on [a.b ], and let G{x) = sSeX 

X ~~ L 

Then 


b n 

D(h) = / G(x)dx -fi£" G(x } ) 


i = o 


= %'(0 + 1? (§5rl G(2M-1) ( a ) - g (2m_1) (m]>* Zm 


( 2 . 12 ) 


+ 0 



as h-*0. 


Proof: We can assume without loss of generality that t— a < b—t m Then, since t is 
one of the Xj , so is b'=2£-a. Furthermore, t is the midpoint o^ the interval [a.b’]. 
Now 


and 


b b 1 


b 


J G(x)dx = J G(x)dx + J G(x)dx 

a a b 1 


(2.13) 


)dr=/ g< x >3< t > -Jc. (2.14) 

O x r 

the integral on the right hand side of (2.14) being an ordinary integral, in which 
the integrand becomes g*(t) when x = £. Applying Theorem 2.1 to the right hand 
side of (2.14), we have 
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Di(h) - *J~ G(x)dx — h £" ~hg'{t) 

Xj L 


zj^b* 

Zj 


= Y 


B 2n 

d 2 ^ -1 

g(x)-g(t) 

(2 At)! | 

dr 2M_1 

x—t 


d 2 *- 1 

g(x)-g(t) 

dx 2 ^ -1 

x—t 


I x=t>' 


A 2 ^ 


+ 0 


^ 2m ] 


as A-*D. 


Now 


and 


V'» 1 

x^ib' Xj-t 

Zj7*t 


■= 0 


d r 

i 

(— l) r 7"! 

dx r 

1 

H 

+ 

-ka 

1 

l 


Combining (2.16) and (2.17) in (2.15), we have 


D\(h) = J G(x)dx — h 2" G{x 5 )-hg'{t) 


Zj^b' 

Zj*t 


77 ^ 1 B s 


_ xr\ -^2a 
m =i (2m)« 


G (2M_1) (a)-G (2 ' u_1) (6') A 2 * 1 + 0^. 2m ] as/i->0. 


Applying Theorem 2.1 to J G(x)dx , we obtain 


D z {h) = fG(x)dx -h 2” G(xj) 


Zj^b' 


= (^ 2 )i G^ 2fx 1 ^(fa') — G^ 2 * 1 ^(b) h z ^ + 0 h Zm as A->0 


M =i (2M)! 

Adding (2.18) and (2.19), (2.12) now follows. 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


[] 
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Corollary: The remainder term 0(h 2m ) in (2.12) is actually given by 

ftaO-.O..*)] = C^\x)dr. (2.20) 

this integral being interpreted as a Cauchy principal value integral. 


Proof: The remainder term 0(/i 2m ) in (2.15) is 


B Zm [(x-a)/h]-Bz m d 2m 

g(x)-q{t) 

(2m)! dx 2m 

x-t 


dx. 


Now making the change of variable x-t +£, we have 


J Bz rn[(x-a)/h]-B Zm d Zm 


Recall that 


(2m)! 


dx Zm .. 


x-t 


dx = 


-L 


&2 m[^/^ l ‘t(^ a )/^] m 

d 2m 

Y 

(2m)! 

df m 



and 


Bk(N+x) = B k (x), N integer, k 2 , 


( 2 . 21 ) 


( 2 . 22 ) 


(2.23) 


B k {-x) = (-1 ) k B h (x) 9 k>2. (2.24) 

Since (t — a)/h is an integer, it follows using (2.23) and (2.24) that the integrand of 

the integral on the right hand side of (2.22) is odd. Consequently, when taken as a 

Cauchy principal value integral, this integral is zero. The result now follows by 

using this in (2.21), and adding to R% m the remainder term of the Euler-Maclaurin 

b 

expansion for the integral J G{x)dx. 

6 * 


[] 


Theorem 2.5: Let g(x) be 2m times differentiable on [a,b], and let 

G(a:) = \x-t \ s g(x). Then 



(2.25) 


= "S' ^-[c»-»(o)-C«»->l<6 )]**» 

- 2 m ^ X - g^\t)h 2 ^ 1 + oM as A-»0. 

^=o 1 * 


Proof: Applying Theorem 2.2 to J G{x)dx, we have 

t 

b 

D x (h) = / G(x)dx -h 2 r G(x y ) 

* xy>< 

- Z> U 1 g M (t )h^ +s+1 + o(/i 2m ) as h ->0. 

, t-a t 

Next applying Theorem 2.2 to J £ s g{t—£)d£ (= J G{x)dx), we have 

0 a 

t 

D z (h) = / G{x)dx -h ^GiXj) 

a xj<t 

= %mr^' Ka)h ^ 

- ^ \-l) M g M(t)h^ s+l + o(/i 2m ) as A->0. 
Adding (2.26) and (2.27), (2.25) follows. 


(2.26) 


(2.27) 


□ 


Theorem 2.6: Let G(x) = \x — t | s log|:n — t \g(x) in the statement of Theorem 2.5, 
everything else being the same. Then 
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D(h) = fG{x)dx -h G{Xj) 


i = o 

Zj*t 


= II§rK ,> <° ) - G ‘ 2 '- ,)( > ) 


h 2 » 


(2.28) 


-2 m ty<'(-s-2ii)+a-^-2^)log h\ 3^£L h 2^l + o^mj as h-*0 


Proof: Similar to that of Theorem 2.5. 


[] 


Corollary: For s=0, (2.28) becomes 

D{h) = g{t)h log h + m S 1 7§frl^- 1 )(a)-G^- 1 )(6)]/i^ 

M_1 ’ (2.29) 

+ 2 ^ 1 ^fe-^r^ (2M) (0 /l2M< ' 1 + 0^ Zm ) as h-> 0. 

Proof: The proof follows by setting s=0 in (2.28) and using the facts that 

£(0) = -1/2, t(“2/2) = 0. fM = 1,2,... , (2.30) 

(see Abramowitz and Stegun (1964, p. 807)). 

[] 

The results in the following theorem will be the ones on which our quadrature 
methods will be based. 


Theorem 2.7: Assume that the functions g{x) and g(x) are 2m times differentiable 
on [tz,b]. Assume also that the functions G(x) are periodic with period T-b —a, 
and that they are 2m times differentiable on H = (— m.oo)\[£ +A:7 T Jjt* = _ 0O . Then 

a) if G(x) = ?(^). and 


<?„[<?] = h t G(xj) . 

J = 1 


(2.31a) 


then 
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E n [G] = [sr(f)+gr'(d)]A + o(/i 2m ) as /i-O. (2.32a) 

b) if G(x) = \x-t |*^(a:) + gr(x), s > — 1, and 


&[G] = h £ G{xj) + ~g{t)h-Zl;(-s)g{t)h s +\ 


i=\ 

Xj^t 


then 


£n[G] = -zY? ^ /p g {Zp) {t )h 2 ^ s+1 + o^ 2m ) 


M=i (2/i)! 

c) if G(x) = |x-f | s log|x-f Iflr(x) + flr(x), s > -1, and 


as /i-»0, 


QnlG] = G(xj)+g(t)h + 2[r(-s)-<-(-s)log/ij ff (f)/i s+1 . 


(2.31b) 


(2.32b) 


(2.31c) 


then 


EnlG] = 2^ 1 [c , (-s-2/x)-^(-s-2 M )log/ij2g|l*2M^^ + 0 |/i 2m ] as h -> 0.( 2 .32c) 
c’) When s=0 in c), by (2.30) and f*(0) = “~~log(27r), (see Abramowitz and Stegun 


(1964, p. 807)), (2.31c) and (2.32c) reduce to 


Qn[G] - A 2 G(x ; ) +g(t)h + log\—\g(t)h. 


£ft 




(2.31c’) 


and 


^ n [G] = 2 m 2 1 ^|^ (z ^0^ 2M+1 + 0^ 2 "‘) as (2.32c’) 


M=1 


where ^[G] = f G(x)dx—Q n [G] in all cases. 


Remarks: 

1) As is seen from (2.32b), (2.32c), and (2.32c*), ^[G] depends only on t , and is 
independent of a and 5. This is a consequence of the periodicity of G(x), and is an 
important property that shall be exploited in the derivation of Romberg -type qua- 
drature formulas in the next section. 

2) Until now we assumed that t is one of the points Xj. When t is arbitrary the 
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periodicity of C(i) can be used to shift the interval [a,b] to [a\b'] such that t 
coincides with one of the Xj in the new interval [a’.b']. This, combined with the 
observation of Remark 1, means that a) E n [G ] in all parts of Theorem 2.7 stays 

n n - 1 

the same if the sum 2] G(xj) in ^.Sla.b.c.c’) is replaced by 2 j G(t+jh) or by the 


identical sum 



G(t +jh), and b) 


Theorem 2.7 holds for all positive integers 


n. These facts will be repeatedly used in Section 4 without further explanation. 


3) In each of the cases of Theorem 2.7, the numerical quadrature formula Q n [G ] is 
computed by using the function values only. The quadrature formula (2.31a) has 
an error which is of order h , and it would seem that one would have to know 
with a high accuracy in order to improve this formula. But as we shall see in the 
next section, the term [p(f )]/i that appears in E n [G ] is easily removed by 
one extrapolation. 

M ^ 

4) Let G(x) = 2fffc(^)log|x"f | k (log\z-t |) p * + 9 \{x)/ (x-t ) + g 2 (x), where a h and 

l 

p k are as described following (1.3), and gk(z ), l<A:<Sf # gi(x) t and gz(x) are 2 m 
times differentiable on [a, 6] and G(x) is periodic with period T=b—a and is 2 m 
times differentiable on Inspection of Theorems 2.4-2. 7 reveals that if we form 
5k [G] as the sum of the quadrature formulas for each one of the terms in G(:c), 
then the error i? n [G] does not contain any contribution from G(x) and its deriva- 
tives at the end points, and the only contribution to E n [G] comes from G(x) and 
its derivatives at x-t , as in Theorem 2.7. 


3. ROMBERG-TYPE NUMERICAL QUADRATURE FORMULAS 

Using the results of Theorem 2.7, we can apply extrapolation techniques to 

b 

derive Romberg-type numerical quadrature formulas for J G(x)dx. 
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The simplest case is that of Theorem 2.7a, and we deal with it first. 

Theorem 3.1: Let G(x) and <? n [G] be as in Theorem 2.7a. Let h k = T/k and 

$n[G] = 2 Q Zn [G] - &[£?] = K f) Gia+jhn-h*/ 2). (3.1) 

3 = 1 

Le., ^[G] is a mid-point rule approximation. Then 

E n {G ] = fG(x)dx - Q n [G ] = okH as ^->0. (3.2) 

a ^ * 

Proof: (3.2) follows directly from Theorem 2.7a, the h term in E n [G] being elim- 
inated when is formed. 

[] 

As a result of Theorem 3.1 we conclude that if G(x) is infinitely differentiable 
on /?=(—' $£*-_«,, then E n [G] tends to zero more quickly than any inverse 
power of n, as We can improve this result considerably whenever G(z) is 

analytic in a strip in the complex z -plane, which, with the exception of the points 
t+kT, A:=0, ±1,..., contains the real line Im z=t/=0 in its interior. 


Theorem 3.2: Let G{z) in the previous theorem be analytic in the strip |Im z | <a, 
except at the simple poles t +kT, A:=0, ± 1 , ±2 Then 




where 


M{t) = max] max | G g (x+ir) | . max | G g (x—ir) | l, 

1— oo<X<w» — »<X<w I 

and 

GAO = |-[g('+0 + G(t-o] . 

Proof: By the periodicity of G{x). we can write 


(3.4) 

(3.5) 
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b t±T/2 

f G{x)dx = f G(z)dx 

a t-T/Z 

77 2 . 

= / 

-r/2 

G 8 (z) is analytic in the strip |Im z | <cr and is periodic with period T . After some 

algebra it can be shown that the n(odd)-point trapezoidal rule approximation or 

772 

the n(ei;en)-point mid-point rule approximation to f G e (£)d£ is just $ n [G]. 

-r/2 

Applying now a theorem due to Davis (1955) (see also Davis and Rabinovitz (1984, 
pp. 314-316)), (3.3)-(3.5) follow. (Actually Davis’ result is stated for the tra- 
pezoidal rule. However, inspection of his proof shows it to be valid for the off-set 
trapezoidal rule. The mid-point rule is an off-set trapezoidal rule.) 

[] 

For arbitrary t , by Remark 2 following Theorem 2.7, the approximation $ n [G] 
can be replaced by 

? n [G]=/>nJ G{t+iK-K/ 2) , (3.7) 

i-i 

and again Theorems 3.1 and 3.2 apply. 

$ n [G] i n Theorem 3.1 has been obtained by employing the Richardson extra- 
polation process once, eliminating the term 0(A) in the error E n \_G\ Since 
j^fG] = 0(A 2rn ) with 771 as large as we wish, there is no need for further extrapola- 
tion. For 9 n [G] as given in Theorem 2.7 b.c.c’, however, we apply the Richardson 
extrapolation process (or generalizations of it) repeatedly in order to eliminate 
successive terms in E n [G\ thus obtaining Romberg-type numerical quadrature 
formulas with higher degrees of accuracy. (A summary of those extrapolation 
methods relevant to the present work is given in the appendix.) 

b 

Let G(x) and £k[G] be as in Theorem 2.7 b,c,c\ and define A = f G{x)dx and 

a 

A (A) = Qn[G\ in the notation of the appendix. Select a sequence of integers 
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\ n i ir=o. 1 ^ n o <ti 1 < • • • , and set h t = r/n { , i =0, 1 Obviously lim h L =0, as 

l~¥ CO 

required; see the appendix. Now the Romberg-type quadrature formulas A$ m \ 
based on/4(hj), 7n<l<7n+q , are of the form 

Af m) = ictfpA(h m+t ). (3.8) 

k =0 

where \ 0 <Jc^q, are constants determined by the nature of the asymptotic 
expansion of A—A(h) = as A-> 0, as described in the appendix. Some of the 

details for two special cases are given below: 

1) If C(x) is as in theorem 2.7b {or Theorem 2.7c’) f then A(h) is of the form given 
in a) of the appendix with r =2, 7 i =s + l+2i, and ft = -2£(-s — 2i)g^ Zi \t )/ (2i)! 

(or r = 2. 7 i = l+2i, and ^i=2^(-2i)g i2i \t)/ (2i)\), i=l,2 Hence for a 

sequence of the form ft = h^p 1 , Z=0,1,... , the can be computed from (A. 6). 
For arbitrary ft, the algorithm given in (A.ll) and (A. 12) in b) of the appendix 
is appropriate with 9 ?(/i)=/l s+3 and r =2 (or <p(h)=h 3 andr = 2). 

2) If G(x) is as in Theorem 2.7c with s?*Q f then A(h) is of the form given in (A.l) 

with e* (h) = [f(“S -2i)— f(— s— 2i)log /i]A 2l+s+1 , and 

ft = 2g iZi) (t)/(2i)\, i = 1.2 The d£) then can be obtained by solving the 

equations in (A. 5). 


Before closing this section we recall that, for any positive integer 
7t, A(/i) = Q n \.G] i n Theorem 2.7 b.^c* has the form 


A(h) = h 2 G(xj) + C(t,h). 

a <i 


;V0 
d +jhtzb 


(3.9) 


where 


C(t,h) = 


g(t)h-2^(-s)g(t)h s+1 

g(t)h+Zl?(-s)-S(-s) log h]g(t)h s + l 


g(t )h 4*log 


2n 


g(t)h 


for Theorem 2.7b 

for Theorem 2.7c 
for Theorem 2.7c\ 


(3.10) 
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4. THE QUADRATURE METHODS FOR INTEGRAL EQUATIONS 

In what follows we consider the integral equation (1.1), with /(£) and g(t) 
being periodic in t with period T , and K{t t x) being periodic both in t and x with 
period T. Of course, a similar treatmet can be given to the eigenvalue problem in 
( 1 . 2 ). 


4.1 The Singular Case 

Let K(t ,x) = Hx{t ,x)/ (t—x) + H z (t ,x). For a given integer N, let 
h = hzN = (6— a)/{2N), and Xj = a+jh, j =1,...,2N . Then setting f=x* for some i, 

b 

and approximating the integral f K(x it x)f (x)dx by the rule ^ in (3.1), we write 

a 

down the following set of equations for the 2 N unknowns Jj (the approximations to 
the corresponding f (xj)): 


where 


uji + 2/ i £ t^K{x i ,x i )Jj = g(xi), i=l,2 2 N, 

i=i 


(4.1) 


s a 


|l if \i—j | odd 
lo if \i—j | even. 


(4.2) 


4.2 The Weakly Singular Case 

We mentioned in the previous section that when G(x) is a known function any 

set of integers 1<71Q<71X< - * • , can be chosen for computing the approxi- 

6 

mation to the integrals J G(x)dx in Theorem 2.7b,c,c\ If, however, we want 

a 

to use the Romberg-type formula for solving integral equations by quadra- 

ture methods, the n t cannot be arbitrary. In fact, we should choose the n* (hence 
the Aj = 7Vn t ) f such that the sets of abscissas that enter the computation of 
-4(^™ +fc) = 0<A:<g-l, where G(x) = K(t ,x)f (x), are all subsets of the set 

of abscissas that enter the computation of A{h m+q ). This is achieved by picking 
the n lt Tn<l^Tn + q — 1, as divisors of ti^^. With this choice of the n t let 
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x j ~ jhm+q - a + l^'^nm+g. With the help of (3.9) it can be verified that, 

for t =Xi, 


n m+7 

= hm+k 2 slj” q ,k G(xj ) 4 * C{zi,h m +k), 

ft 


(4.3) 


where 


1 if \i-j\ is divisible by 
0 otherwise. 

Thus (3.8) becomes, for f =x if 


e m.q.k = 


(4.4) 


. . n m+7 

4 m) = E 


j.=i l*=o 




G(^) + 2 dqWcixi.hn+t) . 


k=0 


(4.5) 


Nowfor^? 4 !, G(xj) = K(x i% Xj)f {Xj). By Theorem 2.7b, c.c’ we note that when 
b) K(t,x) = H j(f ,x ) | £ x | s + H z (t,x), g(x) = /f 1 (f,x)/(x) andg(x) = H z (t,x)f(x) in 
(2.31b). Thus C(f./i) = C(t,h)f(t), where 


C(£,A) = h 


H z (t,t)-Zt(-s)H 1 (t,t)h* 


(4.6b) 


c) K(t ,x) = i/j(f ,x) | f —x | s log 1 1 -x | + H z (t ,x), p (x) and p(x) in (2.31c) are as in b) 
above. Thus C{t,h) = C(f,/i)/(f), where 


(4.6c) 


C(t,h) = h\H z {t,t) + 2[^(-s)-^(-s)log h]H 1 (t ,t)h s 


c’) K{t ,x ) = Hi(t,x )log \t—x \ + H z {t,x), g(x) and g(x) in (2.31c’) are as in b) and 
c) above. Thus C{t,h) = C(f.h)/(f), where 


C(t,h) = h 


H z (t.t ) + log 


277 


Hxit.t) 


(4.6c’) 


Combining the above in (4.5), approximating the integral J K(t ,x)f (x)dx in 

a 

(1.1) by A$ m \ and replacing the f (%j) by the corresponding approximations 
J we obtain the appropriate quadrature methods for (l.l), which are 
defined by the systems of linear equations 
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l<i<n. 


i- 1 


m+qr » 


where 




2 

lfc=0 


Kixi.Xj), j*i , 


(4.7) 


(4.8) 


and 


%l = £ c^W^®,.****). (4-9) 

Jt=0 

with C(£,A) as defined in (^Sb.c.c*). 

It is not the purpose of this work to give precise error bounds or convergence 
results for A = max|/ (xj )— fj | . However, in general, we would expect A to be of 

the order of magnitude of the error in the numerical quadrature formula used in 

b 

approximating the integral J K{t ,x)f (x)dx . Thus, for the singular case, if K(t % z) 

a 

is meromorphic in the strip \Im z | <a with its only poles at f +kT, k = 0, ±1, ±2,... , 
and g(z) is analytic in the same strip, we would expect /(z) to be analytic in this 

strip too. Therefore, by Theorem 3.2, f G{x)dx — Qn[G] = o[e” 27rArff/7T ] as //-»», thus 

a 

we would expect A= o[e“ 27T ^ ff/7T ] as too. For the weakly singular case, by the 

b 

appendix, J G(x)dx — = o|e ?+1 (/i m )j as in general. Thus, we would 

expect A= o|/i^ +1+2? (log as 77t-»°o, where s = max^ a k and p = maxj3 fcl (cf. 

(1.3)). 

Finally, we could use the approximations Ji to /(x*), i= 1 N, where N = 2 N, 

for singular equations and N = for weakly singular equations, to construct a 

trigonometric interpolation polynomial in cos(27rA::r/ T), 

sin(2rrA:x/ T) t k =0,1,... , satisfying Pfii^i) - ^ = 1 A, thus obtaining an approxi- 

mation to / (x) for all x in [a, 6], 
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5. NUMERICAL EXAMPLES 


5.1 The Singular Case 

Example 5.1.1 (See Mikhlin (1964, pp. 122-124)) 

2tt 


o /(0 + 


Z—t 


f (x)dz = it(f ). 


(5.1) 


When a?*0 and a 2 +b 2 ^ 0 (a and b may be complex) a unique solution exists and is 
given by 


nt) = ^TPr“ (t> - 2tt(a*+i,*) j" ( * )<!0t 
b 2 2 r 


z—t 


27ra( 


We first observe that the kernel function K(t,z) = — — cot 

27T 


z—t 


(5.2) 


is periodic with 


period 2rr in both z and t. Also for fixed t , jf(f ,z) is meromorphic in the whole z- 
plane with simple poles at t +2nk , A: =0,±l,:t2,... , thus being of the form described 
in Section 4.1. Next, we observe that if u(t ) is periodic with period 27 r, then so is 
/(f). Also if ix(z) is analytic in a strip | Im z [ <cr, then so is / (z). The last two 
assertions can be verified with the help of (5.2). 


In our numerical experiments we chose u{t ) = (D + cos f)” 1 , D> 1, so that both 
u(z) and /(z) are analytic in the strip |Im z | < a - lo g(D + V D 2 —l ). With this 
choice of Tz(f), (5.2) becomes 


/(0 s 


1 

1 

b sin f 

, fe 2 

a 2 +b 2 1 

D+cos t 

{ Vi? 2 - ij 

a Vi? 2 — l 


. A /V 

Denote 27V, the number of abscissas in (4.1), by N. Then h = 27T/ N and 
Zj = jA, 1 <j < TV. Denote /ftj = / J f and let Aft = ma^|/ (*/)-/ ftj | • 

Then, by what has been said in the paragraph following (4.9), we would expect to 
have Aft = 0(e~ CT ^ /2 ) as N -* This is born out by the numerical results, where, for 
each N, cr can be estimated by the formula 
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a sa 


N m „-N 


•log 


Aft 

A ft, 


- s w, 


'max 


(5.3) 


'max 


where / / m „ is the maximum of the ft s used in the computation. This estimate, of 


course, is based on the Aft, which in turn are based on / (x). An estimate based 
solely on the computed values Tft.; can be obtained from the formula 


a « jlog 


/ ft+2/ ,ft+2 l ~7ft+J.ft+/ 

7KwJ)W~7k.K 


= ffft 


This follows from the following expected behavior of f K.K ' 


(5.4) 


7ft,ft~ f ( x ft) + C( x ft) e as ft -> <» . (5.5) 

rvj .n. no 

Here, of course, /ft,ft and xn=2tt can be replaced by /ftj(ft) and x respectively, 
where Xj(K) = x is the same for all ft used in the computations (in (5.4) they are 
ft, ft+J, and ft+2J). 


Table 5.1.1 gives the results obtained for A#, oft t ft mta * and aft, 
N - 4(4)44, = 44, with a-b-1. Note that 0-44,44, cr 4 and a iz are not defined. 

Note also that, for D= 2, &ft.ft mAX and aft deteriorate for N large. This is due to the 
fact that there is a loss of significance in the arguments of log in (5.3) and (5.4), 
which is caused by the high accuracy of the Jftj. 


A T 


D=l.l 



D=2 


TV 

Aft 

ft 

iV * JY m ax 

on 

Aft 

aft ft 

oft 

4 

2.03x10° 

0.426 


6.10X10 -2 

1.3136 


8 

1.12X10° 

0.4407 


4.60X10- 3 

1.3157 


12 

4.93X10" 1 

0.0445 

0.67 

3.37X10 -4 

1.3170 

1.354 

16 

2.01X10 -1 

0.4442 

0.52 

2.41X10 -5 

1.31677 

1.3195 

20 

7.98X10 -2 

0.4411 

0.474 

1.73x10"° 

1.31683 

1.3171 

24 

3.33X10 -2 

0.4420 

0.456 

1.25X10 -7 

1.3171 

1.316971 

28 

1.39X10" 2 

0.44339 

0.4485 

8.94X10 -9 

1.3169586 

1.3169588 

32 

5.73X10 -3 

0.44303 

0.4456 

6.42X10 -10 

1.316996 

1.3169580 

36 

2.33X10 -3 

0.4400 

0.4444 

4.62X10 -11 

1.3173 

1.3169587 

40 

9.72X10" 4 

0.4420 

0.44391 

3.31X10 -12 

1.3172 

1.316948 

44 

4.01X1Q- 4 


0.44371 

2.38X10 -13 


1.3176 


Table 5.1.1 : Results for and aft, A r =4(4)44, A r mai =44, 

for Example 5.1.1. Exact values of a are 0.44357 for D=l.l and 1.3169579 for D=2. 
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5.2 The Weakly angular Case 

In both of the examples below the kernel function is of the form described in 
c’) of Section 4.2, namely K(t,x) = // 1 (f,z)log|f-z| + Hz(t,x). Therefore, the 
approximations J j to the solution to the integral equation, satisfy (4.7)- 

(4.9) with (4.4) and (4.6c’). Here the O^Ar^g, in (4.8) and (4.9) can be deter- 
mined from (A.5) with e t (A) = A 81 * 1 , i = 1,2 If we let n t =2 l , thus 

h t = T/2 1 , Z = 0,l,2 as we do in the examples below, then the d£ $ are indepen- 
dent of T7i, and can be computed recursively from (A.6) with a \ = 2 -ai_1 , i=l,2 

Thus we find = 1, dfg) = -1/7, = 8/7, dj>$ = 1/217, djfp = -40/217, 

d|[§) = 256/217, etc. From what has been said in paragraph following (4.9) and 
from (A.7), we would expect to have = 0(ffj l +1 ) = 0(2 _ ( 2?+3 ) m ) as m->«, where 
Ai m ^ s max |/ (xj)-fj \ . This is indeed observed numerically. 


Example 5.2.1 (See Christiansen (1971)) 


2tt 


/log 


2 a sin 


ILzeJL 


/ ( x)dx = -^-cos 2f. 


Provided a 5 * 1 , the unique solution to this equation is /(f) = cos 2 1. Otherwise, the 
solution is f (t) = cos 2£+c, c being an arbitrary constant. We observe that the 

[ t —x | 


kernel K{t ,x) - log 


2a sin 


is periodic with period 277 in both t and x, and 


is of the form described in c’) of Section 4.2, namely, 
K(t,x) = Hi(t,x)\og \ t-x | + H z (t 9 x) t with = 1 and H z (t ,t) = log a. 


Table 5.2.1 gives some of the results obtained with a=VT for 
N = n m + q - 2 m+q , m +q = 3(l)7. Note that, for a given row in this table, N t 
the number of abscissas in (4.7)-(4.9), is he same for every member of this row. 
Thus the first column is the result of no extrapolation, the second, of one extrapo- 
lation, the third, of two extrapolations, etc. As can be seen for a given number of 
abscissas N=2 r , the best results are obtained roughly for q &r/ 2. 
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m+q 

0 

1 

2 

3 

4 

5 

6 

7 

3 

3.8X10- 8 

9.9xi n ~ 3 


4.9X10" 3 





4 

4.7X10 -3 

2.3X10" 4 


3.7X10" 4 

4.7xl0~* 




5 

5.0X1O -4 

6.9X10" 8 

4.3x1 o -7 

1 ax If)" 7 

8.BX10" 7 

1.1x10"® 



6 

7.3X10 -5 

2.1X10 -7 

3.2X10" 9 

2.1X10" 10 

fl Qv 10-11 


6.3xlO" 10 


7 

9.2X10 -8 

6.6x10"® 

2.5X10" 11 

5.0X10" 13 

1.5X10" 13 

1 Rx 10-13 

1.8X10" 13 

1.8X10" 13 


Table 5.2.1 : Results for for Example 5.2.1. 

Here N= 2 m * q is the number of abscissas in the quadrature method and q 
is the number of extrapolations in the corresponding numerical quadrature 
formula. The best result for fixed N (or m+q) is underlined. 

Example 5.2.2 (see Henrici (1979, pp. 492-494)) 

Let r:z=z(r), 0<r^ t be the boundary curve of a Jordan region D containing 
the point z=0, and let /(z) be the function mapping D conformally onto the unit 
disk \u\<l in such a manner that /(0)=0 and /'(0)>0. To determine /(z) it is 
sufficient to know its values on the boundary of D . Because |/ (z) | =1, it suffices to 
know t5(t) = arg(/ (z(r))), 0^r<£. Let £(r) = ^'(r). Then, provided the capacity of V 
is different than 1, £(t) is the unique solution of 
P 

f log I z (ct)-z (r) I £(r)dr = 2 tt log | z (a) \ , 0<<r</3, (5.6) 

o 

which is known as Syram’s equation (see Symm (1966)). The capacity of V is 
different than 1, in particular, if T is entirely within or entirely without the unit 
circle. 

Obviously, the kernel AT(<7 ,t) = log | z (ff)-z (r) | is periodic in both a and r with 
period and is of the form K(a t r) = ffi(cr,T)log | a—r\ 4- Hz(<J,t) with Hi(o,a) - 1 and 
= log | z'(a) | . Also, log | z (r) | and £(r) are periodic with period /?. 

Let /?=2 tt and let D be the elliptic domain whose boundary curve is 
V: z(t) = c(e lT + se _lT ) l (Kt<27T i with 0<s<l, and c >0 chosen so that V is entirely 
without the unit circle. (Actually the capacity of T in this case is c , so that it is 
sufficient to choose c / 1.) The semi-axes of D are c(l+£) and c( 1—e). The solu- 
tion for f(r) can be expressed as 
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£(r) = 1+4 £ (-i) fc ( 2fcT )- 

Note that both log|z(T)| and £(t) are analytic functions of T. 

Observe that £(r) for this example is symmetric with respect to both the Re z 
and the Im z axes. This can be utilized to reduce the dimensions of the matrix by 
4, thus reducing the storage and computing time considerably. 

Tables 5.2,2a and 5.2,2b give some of the results obtained with c=50 and 
g=0.1 and £=0.5 respectively, for N =7^^ = 2 m+q , m+q <7. As in Table 

5.2.1, in these tables too, for a given row, N , the number of abscissas in (4.7)- 
(4,9), is the same for every member of this row. 


77L + q 



5 . 

0 

1 

2 

3 

2 

1.6X10" 1 




3 

2.9xl(T z 

2.7X10 -2 



4 

4.0X10' 3 

wmmsm 

4.5X10 -3 


5 

5-OxlO- 4 

2.7X10" 5 

6.1X10- 5 


6 

6.3X1CT 5 

7. lxlO -7 

l.OxlO -7 

4.8X10" 7 

7 

7.8X10" 6 

2.2X10 -8 

6.7X10 -10 

1.5X10 -10 


Table 5.2.2a : Results for for Example 5.2.2 with £=0.1 and c =50. 

Here N= 2 m+q is the number of abscissas in the quadrature method 

and q is the number of the extrapolations in the corresponding numerical 

quadrature formula. The best result for fixed N (or m+q) is underlined. 


77L~bg 


9 


0 

1 

2 

5 

5.7X10" 2 

1.9X10 -1 


6 

1.6X10" 3 

1.5X10 -3 


7 

9.4x 10 -4 

3.2X10 -5 

3.6X10 -5 


Table 5.2.2b : Results for for Example 5.2.2 with £ = 0.5 and c =50. 
Here N=2 m + q is the number of abscissas in the quadrature method 
and q is the number of extrapolations in the corresponding numerical 
quadrature formula. The best result for fixed N (or m+g) is underlined. 


For small values of £ the ellipse is close to a circle. Therefore, £(r) does not 
change rapidly with a, and this explains the high accuracy obtained for the 
approximations to £(t) even with a small number of abscissas when £=0.1. For 
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large values of £, however, the ellipse is elongated, and this leads to rapid changes 
in £(t) in the vicinity of r = rr/2 and r= 377/2, Le., where £(r) is maximal. This 
then explains the slow convergence of the approximations for £=0.5. Furthermore, 
extrapolation becomes effective in this case starting with a relatively large N . 

To improve the performance of the quadrature method above for the cases in 
which f(r) has rapid changes we can make a change of variable of integration 
T = r(Y0 so that £( T (V0) changes slowly as a function of This can be achieved by 
picking t(^) such that dr/ d^ becomes small where d£/dr is large. This is 
equivalent to having more abscissas in places where £(r) changes rapidly. Need- 
less to say, the transformation r = t(Y 0 should be such that dr/ dip is a periodic 
function of 


For the example under consideration we can choose 


dj/ __ 77 V I+ 77 2 
dr ?7 2 +cos 2 t 
so that 


rj a positive constant, 


^(t) = tan 


-1 


Vl+ 77 2 


V 


■tan r 


or 


0<t<tt/ 2 


t(^) = tan 


-1 


JL 


-tan m 


(Vi+tj 2 
r(^) is now extended so that 


0 2 . 


r(^) = 7T - t(77 — Tp), 77/ 2<Tp<7T , 
r(^) = 7T + r(^“ 77), t<t/^2tt . 

We used this transformation with different values of 77 . In Table 5.2.3 we give some 
of the results with c =50 and £=0.5 obtained for the errors at r = n/ 2, the point at 
which the error is maximum, and at r= 0 . The number of abscissas in all cases is 
jV= 32 , and g= 0 , i.e., no extrapolation is employed. Nevertheless, the improve- 
ment in the results is remarkable. 
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V 

error for r=0 

error for t=tt/ 2 

0.2 

0.5 

100 

4.0 x 10 -3 
1.2 x 10~ 4 
2.4 x 10 -4 

8. lx 10 -4 
2.3 x 10" 3 
5.7 x 10 -2 


Table 5.2.3 : Results for the error at r=0 and t= 7T/2 in the approximations 
to f(r) in Example 5-2.2, with c =50 and £=0.5, and the change of variable r=r(^) 
as described in the text. The numerical quadrature formula used has 
TV =32 abscissas and does not employ extrapolation. 
f(0) = 0.014671... and ^(tt/ 2) =4.5324.... 


6. CONCLUDING REMARKS 


In this section we shall briefly discuss some known quadrature methods that 
are related to those proposed in the present work. 

A lot of attention has been paid to Cauchy principal value integrals and singu- 
lar integral equations. We do not intend to survey all the methods developed for 
these, but we shall restrict our attention to those that are periodic. When 


W.X) = 


t—x 


, Q<t 1 2<27r 1 (the Hilbert kernel), 


Kij = 2hSijK(x i9 Xj) = i-cot 


(i-j)n 
2 N 


ij=l 2 N, 


(cf. (4.1)). The matrix = (X* ; ) is called Wittich’s matrix (see Gaier (1964, p. 76)) 
and its properties are well known. Gutknecht (1981) has used Wittich’s matrix to 
discretize Theodorsen's integral equation for conformal mapping and has 
analyzed various nonlinear iterative techniques for solving the resulting equa- 
tions. We note that Theodorsen’s equation has the Hilbert kernel as its kernel. 


The Hilbert kernel arises as part of the Cauchy kernel on closed curves and 
Atkinson (1972b) has proposed and analyzed product-type integration formulas 
for the Cauchy transform 
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f r closed, z e T, 

which are different than those proposed in the present work. 

As for the weakly singular integral equations, the case that has received the 
widest attention is that of logarithmic singularity. Periodic integral equations 
with logarithmically singular kernels arise naturally, for example, in conformal 
mapping (Symm’s equation) and two-dimensional potential theory. Two of the qua- 
drature methods that have been considered for such equations are based on the 
so called modified quadrature method (see Kantorovich and Krylov (1984, p. 102), 
and Baker (1977, Chap. 5, Sect. 4)). In this method we begin by writing (l.l) in the 
form 


b b 

af{t) + fK(t.x)[f(x)-f(t)]dx + f (t)f K(t ,x)dx =g(t) . (6.1) 

a a 

Employing the numerical quadrature formula ^ Wj F(Xj ) to approximate the 

i = i 

6 

integral f F{x)dx , and using the fact that lim K(t ,x)[f (x)— / (t )] = 0, we replace 

a 


the integral equation above by the linear equations 


CJ + f K(z il x)dx 


fi + =^( :c i). i = l.— 

j=i 


( 6 . 2 ) 


where Ji are approximations to the /fo). Kussmaul and Werner (1968) have 

applied this method with equidistant abscissas Xi+i—Xi = h, 'l = 1,...,ti—1 # to periodic 

integral equations with logarithmically singular kernels and have shown that the 

b 

error is 0(/i 3 ) as assuming f K(t i3 x)dx has been computed exactly. For ker- 


nels of the type K(p,q) = log p(p,q) t where p(p,q) = ^(x (p )~x (q)) 2 +(y (p)-y(q )) 2 , 

and (x(p), y(p)), a<p<b t is the parametric representation of a simple closed 

curve in the x-y plane, Christiansen (1971) has used the modified quadrature 

b 

method in (5.1) with J K(t it q)dq essentially replaced by a numerical quadrature 

a 
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approximation. The numerical results indicate that this method too has an error 
of 0(/i 3 ) as 7i-»0. For both methods above we can show, using Theorem 2.7c\ that 
the errors in the numerical quadrature formulas are 0(7i 3 ) as A-*0, although this 
proof for Christiansen’s method becomes very complicated. 

For weakly singular Fredholm integral equations of the second kind methods 
based on product integration have also been developed and analyzed by Atkinson 
(1967, 1972a). 

Finally, the approach of this work can easily be extended to coupled 
Fredholm integral equations for several unknown functions in which several types 
of singularities occur simultaneously. 
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APPENDIX 

In this appendix we summarize some of the important aspects of generaliza- 
tions of the classical Richardson extrapolation process. 

Let the function A(h), where A>0 is a continuous or discrete variable, have 
the asymptotic expansion 

A(h) ~ A + f) ) as /i-»0, (A.1) 

i=l 

with e i+1 (/i) = as A-*0, £ = 1,2 A(h) and the e^h) are known, but A and 

the & are not. We are interested in approximating A , which, in general, is limA(A) 

h -*0 

when this limit exists. 

Select a sequence of h’s, namely h 0 >hi>... , such that lim/ij = 0, and define 

l 

AyP (the approximation to A) and the & to be the solution of the system of linear 
, equations 

A(h t )=A^U 3—l—j +n - (A.2) 

1 = 1 

For some classes of functions A(h) and e i (/i),i=l,2 and some choices of 

h L , Z = 0,1,... , the following can be shown: 

1) For fixed n, 

A-AV) = o|e 7l+1 (/i J -)j as j (i.e. , hj-*0), (A.3) 

2) For fixed j , A$^-*A as n-»oo f and the convergence in this case is better than in 

the previous case. 

For the present work, it is important to note that can be expressed as 

EdV}A(h j+h ), (A.4) 

Jfe=0 

where the coefficients can be obtained directly by solving the linear system of 
equations 
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£d&=l 

k n „ (A.5) 

k= 0 

The problem of computing the recursively has been attacked by several 
authors, Schneider (1975), Havie (1979), and Brezinski (1980) have devised an 
algorithm that has been denoted the E-algorithm. Recently, a more efficient algo- 
rithm has been derived by Ford and Sidi (1984). 


Special Cases 

a) e i {h) = h 7 \ 0<yi<y 2 < .... limy* = 

i -* « 

For the choice h t = /i 0 p*, £=0,1,... , 0<p<l, the can be computed from the 

recursion relation 


<03 = 




(T n ~l 


0<k<n, 


(A- 6) 


with = 0 = c^'^ +lf 7i=0,l =1, ;=0,1,.„ , and <Ji=p 7i , £ = 1,2 Obvi- 

ously, the are independent of j. This development, in slightly different nota- 


tion, is due to Bulirsch and Stoer (1964), who also give a recursive algorithm for 
and a thorough convergence analysis. First, 


A~A = 0(ctA+i) as 

Second, if there exist constants &=1,2,... , for which 


| A(h)-A-'£p i h 7i \4~fyh N , h<h 0 , 

i = 1 


and, for some fixed i?>0. 


(A-7) 


(A-B) 


& =0 (k\R k ) as (A.9) 

then, for the special case yi = yo+ir, for some 7o>0 andr>0, 

A— 4p) = 0(cj nB ) as n->°o, (A- 10) 

where o=p r/2 + £<1, any £>0. For a similar result pertaining to the case of arbi- 
trary y i see Bulirsch and Stoer (1964). 
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b) e t (/i) = <p{h)h*-r % r> 0. 

For arbitrary h it the can be computed by using the recursive W-algorithm of 
Sidi (1982). As a consequence of this algorithm we have 



i=0 


□<^< 71 , 


where can be obtained from the recursion relation 


(A- 11) 


m = ^ . CM^n, 


hf+ n -hj 


(A. 12) 


with 6$-i = 0 = n=0,l,... , and <5<$ = l/<p(hj), j = 0,1 We note that this is 

also a special case of the generalized Richardson extrapolation process of Sidi 
(1979), a very efficient recursive algorithm for which has recently been given by 
Ford and Sidi (1984). 


We also note that the algorithms given for the two special cases above are 
more efficient than the algorithms for the general extrapolation algorithms for 
obtaining A ^ defined by (A.2), since they take advantage of the special forms of 
the problem. 
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